Energy resolution and discretization artefacts in the numerical renormalization group 
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We study the limits of the energy resolution that can be achieved in the calculations of spectral functions 
of quantum impurity models using the numerical renormalization group (NRG) technique with interleaving 
(z-averaging). We show that overbroadening errors can be largely eliminated, that higher-moment spectral 
sum rules are satisfied to a good accuracy, and that positions, heights and widths of spectral features are well 
reproduced; the NRG approximates very well the spectral-weight distribution. We find, however, that the dis- 
cretization of the conduction-band continuum nevertheless introduces artefacts. We present a new discretization 
scheme which removes the band-edge discretization artefacts of the conventional approach and significantly 
improves the convergence to the continuum (A — > 1) limit. Sample calculations of spectral functions with high 
energy resolution are presented. We follow in detail the emergence of the Kondo resonance in the Anderson im- 
purity model as the electron-electron repulsion is increased, and the emergence of the phononic side peaks and 
the transition from the spin Kondo effect to the charge Kondo effect in the Anderson-Holstein impurity model as 
the electron-phonon coupling is increased. We also compute the spectral function of the Hubbard model within 
the dynamical mean-field theory (DMFT), confirming the presence of fine structure in the Hubbard bands. 

PACS numbers: 71.27,+a, 05.10.Cc, 72.10.Fk, 72.15.Qm 



I. INTRODUCTION 

Condensed-matter systems often exhibit rather complex be- 
havior due to strong Coulomb repulsion between the elec- 
trons at short distances. These effects become very pro- 
nounced when electrons are strongly confined either in in- 
ner electron shells (transition and rare-earth atoms) or in ar- 
tificial nanostructures (quantum dots). Theoretical studies of 
the corresponding many-particle problems rely increasingly 
on advanced computational techniques such as the numeri- 
cal renormalization group (NRG)^ 2 j 3 -. The NRG allows to 
study both static and dynami c 4 ' 5 ' 6 ' 7 ' 8 ' 9 ' 10 properties of quan- 
tum impurity models like the Kondo model or the Ander- 
son impurity model. Applications range from studies of ther- 
modynamic properties of magnetic impurities in norma l 2,11 ' 12 
and superconducting 13 ' 14 host metals, dissipative two-state 
systems 15 , electron transport through nanostructures 16 , to the 
use of the NRG as an impurity solver in the dynamical mean- 
field theory (DMFT) 17 ' 18 ' 19 ' 20 . 

The foundation of the NRG is the transformation of a model 
with an infinite number of degrees of freedom (the continuum 
of the conduction-band electron states) to a model with a finite 
number of lattice sites (known as the "hopping Hamiltonian" 
or the "Wilson chain") which is numerically tractable using 
a computer. This transformation consists of three steps: 1) 
logarithmic discretization of the conduction band into increas- 
ingly narrow intervals around the Fermi level, 2) dismissal of 
combinations of states which do not couple directly to the im- 
purity, and 3) unitary transformation to a basis in which the 
conduction-band Hamiltonian takes the form of a semi-infinite 
chain with exponentially decreasing hopping between neigh- 
boring sites. In the first step, the discretization is controlled 
by a parameter A > 1, which sets the energy widths ~ A~™ 
of the intervals; the continuum is restored in the A — > 1 limit, 
while typical values used in practical calculations are A = 2 



or even much higher, depending on the application. The main 
approximation in the NRG intervenes in the second step (dis- 
missal of higher modes); this approximation is controlled and 
it becomes better as A is decreased 1 . An alternative discretiza- 
tion scheme 21 leads directly to the decoupling of higher modes 
at the price of using a non-orthogonal basis. The third step 
(mapping from the "star Hamiltonian" to a "chain Hamilto- 
nian") can, in fact, be omitted 22 at the cost of significantly 
higher computational requirements. 

After these initial steps, the Hamiltonian is diagonalized 
iteratively, taking one more chain site into account in each 
NRG iteration. Since the Hilbert space grows exponentially, 
only a finite number of low-lying states are kept in each itera- 
tion, while high-energy states are discarded (truncated). This 
procedure is possible due to the "separation of energy scales" 
which simply means that the matrix elements between the bot- 
tom and top end of the excitation spectrum are small 1 ; this is 
an important property of quantum impurity models. Trunca- 
tion is another source of systematic errors in NRG. These er- 
rors are more difficult to estimate a-priori, but they can be kept 
small by a proper choice of A and by performing the trunca- 
tion at suitably high cutoff energy. 

While the NRG is the method of choice to study low-energy 
properties of quantum impurity models, it is, however, com- 
monly believed that it has inherently limited energy resolution 
at higher energies due to the discretization of the conduction 
band. This is particularly relevant for the calculations of dy- 
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such as the impurity spectral function 



or the dynamical susceptibilities. Since the continuum impu- 
rity model is mapped onto a finite chain, the spectral func- 
tion consists of a set of delta peaks with given energies and 
weights. These peaks need to be broadene d 3 ' 10 ' 23 to obtain 
the desired final result: a smooth spectral density function. In 
order to efficiently smooth out spurious oscillations, broaden- 
ing kernel functions with long tails are usually chosen. The 
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log-Gaussian broadening function exp (—(In a; — lnu/) 2 /6 2 ) 
is very commonly used since it is well adapted to the loga- 
rithmic discretization grid. Unfortunately, the slowly decay- 
ing tails lead to strong overbroadening effects, restraining the 
effective energy resolution at higher energies and completely 
washing out any narrow spectral features with small spectral 
weight. 

Narrower broadening functions can be used when the so- 
called interleaved method (also known as the "z-averaging") 
is use d 21 ' 24 ' 25 ! 26 . The interleaved method consists of perform- 
ing several NRG calculations for different (interleaved) log- 
arithmic discretization meshes controlled by the "twist" pa- 
rameter z G (0 : 1], In this way, the information is sampled 
from different energy regions in each NRG run. The spectral 
function is then computed by averaging over all z values. Al- 
though the interleaved method does not truly restore the con- 
tinuum A — ► 1 limit, it is surprisingly successful in removing 
oscillatory features in the spectra; even averaging over only 
two values of z is often very beneficial. 

In this work, we study to what extent the energy resolu- 
tion of the NRG can be ultimately improved by the interleaved 
method. We perform the averaging over a very large number 
of values of z and use very narrow Gaussian broadening ker- 
nel of width proportional to the energy of each individual delta 
peak. This approach, although rather costly in terms of the 
required computational resources, eliminates overbroadening 
and provides spectral functions with very high energy resolu- 
tion even on the energy scale of the width of the conduction 
band. In addition to allowing us to study the fine structure in 
the spectral functions of impurity models, this high-resolution 
approach also uncovers the artefacts which are inherent in the 
NRG and cannot be eliminated by the z-averaging. The arte- 
facts diminish as A is decreased, but they are present in any 
practical NRG calculation. By determining the appearance of 
the artefacts and their expected locations, one can properly 
take them into account when interpreting the results. We also 
propose a new discretization procedure which is very success- 
ful in removing the most severe NRG discretization artefacts. 
This improvement makes NRG a powerful technique for ac- 
curately studying both low and high energy scales, thereby 
increasing its value as a reliable impurity solver in DMFT 

This work is structured as follows. We introduce the Ander- 
son impurity model in Sec. [II] and the details of the NRG cal- 
culations in SecHIIl To explore how accurately NRG approx- 
imates the spectral-weight distribution, we present in Sec.lIVI 
the sum rules for spectral functions of the Anderson impu- 
rity model, the fulfilment of which is then studied in Sec. [V] 
The discretization artefacts are discussed in Sec. |VI| while 
in Sec. I VIII we present the modification to the discretization 
scheme which renders these artefacts less severe. In Sec. lVIHl 
we present examples of high-resolution spectral functions for 
the Anderson and Anderson-Holstein impurity models which 
reveal interesting details, which cannot be easily obtained by 
any other method. Finally, in Sec.|IX]we demonstrate the fea- 
sibility of using the high-resolution NRG approach in a DMFT 
setup. The resolution is sufficient to resolve the fine structure 
in the Hubbard bands, in particular the accumulation of the 
spectral weight at inner Hubbard band edges. 



II. ANDERSON IMPURITY MODEL 

We consider the Anderson impurity model 27 , the paradigm 
of the quantum impurity models. It is defined by the following 
Hamiltonian 

( 1 ) 

where operators c ka describe the continuum conduction-band 
electrons and operators d c the impurity level, e k is the band 
dispersion, V k the impurity hybridisation, TV the number of 
the lattice sites, e the impurity energy and U the on-site 
electron-electron repulsion. Furthermore, n a = d\d a and 
n = n-f + ni. I n me derivations to follow, it is more con- 
venient to rewrite the hybridisation part of the Hamiltonian 

H hyh = VJ2(fld°+4fo*)- (2) 
Here the hybridisation constant V is defined as 

k 

and the operator /q ct as 

The operator fo a thus describes the combination of band 
states which couple directly to the impurity level. The hy- 
bridisation strength is given by T = irpV 2 , where p is the 
density of states (DOS) in the conduction band. In numerical 
calculations we will use a constant DOS p = 1/2D, where 
2D is the bandwidth, unless noted otherwise. 

In the NRG, the continuum of band electrons is reduced to 
the hopping Hamiltonian 

oc 

C=E^(iw+ H 4 (5) 

i=0,<T 

The operator /o CT represents the previously introduced combi- 
nation of states, while fi„ for i > 1 describe further orbitals 
along the Wilson chain. The coefficients ti depend on the dis- 
cretization scheme and on the parameters A and z\ asymptoti- 
cally they behave as t j ~ A~ J /2 . We emphasize that this is not 
an exact representation of the continuum band Hamiltonian. 

m. METHOD 

Dynamical NRG calculations were performed using the 
density-matrix approach 4 - 28,29 using the density matrix com- 
puted at the energy scale of 10~ 12 D. Spectral functions were 
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obtained by delta-peak broadening using a Gaussian kernel 
with a width proportional to the peak energy^ii: 



P(cj,E) = 



1 



(6) 



where u> is the energy of the point in the spectrum, E is the 
delta-peak energy and the width of the Gaussian is T]e = v\E\ 
with r\ a constant (we mostly use r\ = 0.01 or r\ = 0.015); 
the relative spectral resolution is thus expected to be constant, 
AE/E rj -q. For the purposes of obtaining high-resolution 
spectral functions, it is very important to use Gaussian broad- 
ening rather than, for example, Lorentzian broadening, due 
to the fast decrease to zero of the Gaussian function. We 
also note that the conventional log-Gaussian broadening ker- 
nel exp (— (lnw — lnti/) 2 /6 2 ) becomes equivalent to a sim- 
ple Gaussian kernel for small enough b, aside from a small 
asymmetry of the log-Gaussian function. Furthermore, pa- 
rameters rj and b are related by b = \/2ij in this limit. Nev- 
ertheless, the symmetry of the Gaussian function is beneficial 
for the purposes of this work. For some further comments on 
the spectral function broadening, see Appendix A. 

The discretization was performed using the non- 
orthogonal-basis-set approach of Campo and Oliveira 21 , 
with averaging over N z = 32 or N z — 64 values of the twist 
parameter z, equally distributed in the interval (0 : 1]. We 
note that in order to obtain a smooth spectrum, r\ and N z need 
to be chosen such that i]N z is of order 1 . 

The truncations were performed at an energy cutoff 
-E-cutoff = IOwjv, where ujn <x At n / 2 is the characteristic 
energy scale at the jV-th NRG iteration. When necessary, ad- 
ditional states were retained above this cutoff energy to ensure 
that the truncation was performed within an energy "gap" of 
at least Q.OIlon, so as not to introduce systematic errors which 
may arise by retaining only parts of clusters of nearly degen- 
erate states. Charge conservation and SU(2) spin invariance 
have been explicitly taken into account. 

Spectral functions were obtained by "patching" together 
spectral functions from every second energy shell (the N/N 
2 approach) 23 . The details of the patching approach are impor- 
tant and, if not done properly, the procedure will accentuate 
the discretization artefacts. At every even- N NRG interaction, 
we perform the patching as described in Ref. 23: we merge 
spectral peaks in the energy range [pu>M : pAusx] (unmodi- 
fied) and spectral peaks in the range [pAuiN ■ pA 2 lon] (after 
linear rescaling) with the total spectral density; p is some con- 
stant that we refer to as the "patching parameter". We return to 
the patching procedure in Sec. [VI] where we also comment on 
the relative merits of the patching approach and the complete- 
Fock-space technique^. 



IV. HIGHER-MOMENT SPECTRAL SUM RULES FOR 
THE ANDERSON IMPURITY MODEL 

A simple way of quantifying the distribution of the spectral 
weight is through the moments, defined as 



Mr. 



where A a {uS) = — — ^xa{{d a ; d^))^ is the spectral function. 
A stringent test of the calculated dynamic property (spectral 
function) is to verify that it satisfies the sum rules which relate 
the moments to various static quantities (expectation values). 
The zero-th moment is simply the normalization condition for 
spectral functions 



Mo = 1. 



(8) 



Higher-moment spectral sum rules for the Anderson impurity 
model can be derived as?^ 

» m = ({[d„,H] m ,4}), (9) 

where [A, B] m is the iterated commutator, defined recursively 



[A,B]i = [A,B] = AB-BA 
[A,B] n+1 = [[A,B] n ,B] 



(10) 



while {A, B} = AB + BA is the anticommutator. The first 
moment (mean energy) is simply the Hartree energy of the 
impurity level, 

Mi = e + C/(n_ CT ), (11) 

while the second is 

M 2 = V 2 + e 2 + (U + 2e)U (n_«r) . (12) 
The variance of the spectral function is thus 

k 2 = M2 - /i 2 = V 2 + U 2 (n_ a ) (1 - (n_ CT )), (13) 

i.e. a sum of the hybridisation width V 2 = T/(irp) and the 
interaction-induced width. The third moment is 

+ ^ 3 =e 3 + 2eV 2 + U{3e 2 + 3eU + U 2 + 4V 2 ) 



UV 
toUv(h W a 



(4V + (tf + 2e)(j 



(0) 



(14) 



Here the operator u/, ff is the /o-orbital occupancy nf t<7 — 

t (i) 

/(kr/oo- and the operators h are hopping operators h a — 
d\fi >a + fi a de between the impurity orbital and the site i 
of the Wilson chain. The third central moment is thus 

«3 = M3 - 3/^1^2 + 2fl\ = 

U 3 (2 (n_ CT ) 3 -3(n_ CT ) 2 + (n_ CT )) 
-V 2 (e + U(2(n f ^)-(n_ a ))) 



(15) 



uj m A a (uj)duj. 



(7) 



which simplifies in the non-interacting limit to K3 = — eV 2 . 
The fourth moment is 
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/i 4 =e 4 + 3e 2 V 2 + V 4 + U (4e 3 + 



6e 2 U + 4eU 2 



U 3 + 2{7e + AU)V 2 



UV 



(U + 2e) 2 (hVty + V ((8e + 3U) (n f ,-*) + U (<?_ CT ))] + t 2 V 2 + 2t U(U + 2e) (/i^) 



(16) 



where operator g„ = T + 2(0± + n<jrif a ); here T = 
dtd[/o,t/o,J. + h.c. is the two-particle hopping operator and 

0± = <i|dx/d j./o,t + n - c - i s tne transverse part of the spin- 
exchange operator. In the non-interacting limit, the fourth mo- 
ment simplifies to 



+ (3e 2 + tl)V 2 + V 4 



(17) 



It is important to point out that the expressions for /13 and 
/i 4 depend on the discretization through the coefficient to and 

the operator h_ a (for /j,3 this is the case only for U 7^ 0). They 
are therefore not exact. While it is possible to derive exact ex- 
pressions in terms of Vk, efe and operators d\ck,<j + H.c., they 
are of little practical use. This implies that in the interacting 
case, calculations of [13 and /14 and the fulfilment of the cor- 
responding sum rules must be considered above all as a test 
of the internal consistency of the method and of the extent 
of errors brought about by the NRG truncation ("energetics"). 
Comparison with exact [j>3 and \i± (were they known) would 
inevitably show some discrepancy (in the following, we will 
demonstrate such behavior for ^4 in the non-interacting case). 



V. SPECTRAL WEIGHT DISTRIBUTION AND SUM 
RULES 

A. Non-interacting case 

We first consider the non-interacting (U = 0) resonant- 
level model. The spectral moments are tabulated in Table U 
The spectral function of this model is given exactly as 



with 



A(u) 



A(w) = r 



1. 



-Im 



1 



AM 



1, (1-uD 
1 + - In — 

7T \l+LO D 



(18) 



(19) 



for oj E [-D, D]. This expression for A(uS) is used to com- 
pute the reference values for spectral moments exactly (sec- 

(e) 

ond column, n\ '). The right-hand sides of the sum rules, 
Eq. (0, are computed in the standard way with j3 = 0.75 
(third column, fj,y')i^. The fourth column contains mo- 
ments calculated by summing the suitably weighted delta- 
peak contributions to the spectral function, ^j!f \ and finally 
the fifth column contains moments calculated directly by per- 
forming a numerical integration with a spectral function after 



(s) 

The first three moments calculated as static quantities, fi\ , 
trivially agree with exact values since they are constants, 
while there is a 7 percent discrepancy for the fourth. This 
can be attributed to the discretization errors as described pre- 
viously in Sec. UV] It must be noted, however, that the fourth 
moment of a Lorentzian peak located near the Fermi level 
strongly depends on the details around the band edges and 
contains little information about the spectral distribution in 
the frequency range of interest (i.e. around the peak itself). 

(s) 

More importantly, we find good agreement between fj}- ' and 



and 



broadening, /1 



(>>) 



the moments computed from dynamic quantities, 

fi\ b \ with errors in the few permil range. This internal self- 
consistency of the method implies that the accuracy of the en- 
ergy levels in the range where the contributions to the spec- 
tral function are sampled from is very good. The difference 
between results from a calculation from delta-peak weights, 

Hi , or from broadened spectral function, iif \ is remarkably 
small. This already suggest that the broadening procedure it- 
self does not lead to any appreciable overbroadening. 

To study how the logarithmic discretization affects the 
spectral weight distribution, we plot the spectral function of 
the non-interacting model for a range of values of the dis- 
cretization parameter A, Fig. Q] The peak position, width and 
height are well reproduced; the position to within less than 
one percent even at A = 2, while the height and the half- 
width at half-maximum both deviate by less than 5 percent. 
As expected, the agreement improves as A is decreased, al- 
though not in a uniform manner. It may be noted that some 
spectral weight seems to be missing in the peak (with the sit- 
uation improving as A -> 1). This is indeed the case; the 
missing spectral weight is located in the NRG discretization 
artefacts that are the topic of Sec.lVIl 



B. Interacting case 

We now switch on the interaction and consider an asymmet- 
ric Anderson impurity model in the Kondo regime, U/irT ^> 
1. Exact results for moments are not available in this case, 
but we can compare /i^ and fJ-[ d \ Table HT1 We find a simi- 
lar degree of agreement (few permil) as in the non-interacting 

case. We also observe that the moments [if^ ( an( l /4 ) calcu- 
lated for each value of z separately depend relatively little on 
z. This is somewhat surprising given that unaveraged spec- 
tral functions are extremely oscillatory. It also implies that if 
we are really interested in a quantity which can be expressed 
as an integral of the spectral function multiplied by some rel- 
atively smooth weight function, there is only little benefit in 
performing the z-averaging. 
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Moment 



Exact, (i\ 



(e) 



Static, (i\ 



Dynamic (delta peaks), (i\ 



(d) 



Dynamic (broadened), fi. 



(b) 



(10 
(11 
(12 
(13 
(14 



1 

-0.050000 
0.0056831 
-0.00044331 
0.00110129 



-0.050000 
0.0056831 
-0.00044331 
0.0010225 



0.999442 

-0.049983 

0.0056866 

-0.00044366 

0.0010220 



0.999981 

-0.049999 

0.0056871 

-0.00044389 

0.0010225 



Table I: Moments for the non-interacting impurity model with parameters e/D — —0.05 and F/D = 0.005. NRG parameters are A = 2, r\ 
0.015, iV z =32,p = 2. 



Moment 



Static, [i\ 



00 



Dynamic (delta peaks), \i 



(d) 



Dynamic (broadened), fi { 



((.) 



(10 

Mi 

(12 
(13 
(14 



-0.0123204 
0.00455271 
-0.000138146 
0.0010179 



1.000303 

-0.0123184 

0.00455556 

-0.000138222 

0.00101737 



1.000306 

-0.0123184 

0.00455549 

-0.0001381970 

0.00101748 



Table II: Moments for the asymmetric Anderson model with parameters U/D = 0.07, e/D — —0.05, F/D = 0.005. NRG parameters are 
A = 2, Tj = 0.015, N z = 32 and p = 2. 
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exact 
A=1.6 
A=1.7 
A=1.8 
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1.0032 

1.0177 

1.03225 

1.0553 

1 .0502 

1.0424 



-0.049 

0.05, r/D=0.005, U=0 - 
_r|=0.015, N =32, p=1.75 



0.048 0.052, 



-0.05 



-0.04 



co/D 



-0.03 



Figure 1: (Color online) Spectral function of the non-interacting 
model for a range of discretization parameters A, compared with the 
exact solution, Eq. l !18t . 



We now study the spectral function of the symmetric An- 
derson impurity model shown for a range of discretization pa- 
rameters A in Fig. |2] The spectral functions overlap to a very 
good approximation and there is little systematic overbroad- 
ening. The width of the charge-transfer peak is, as expected, 
approximately 2T. The Kondo resonance is well reproduced 
with a notable exception of A = 1.8, where we find an arte- 
fact which takes the form of a depression at the top of the 
Kondo resonance. For this value of A, the Friedel sum rule 
A(ui = 0) = l/Vr is strongly violated. This is another man- 
ifestation of the NRG artefacts that will be discussed in the 
following; the result is improved by tuning the patching pa- 
rameter p. 

A very successful method to reduce overbroadening ef- 
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Figure 2: (Color online) Spectral function of the symmetric Ander- 
son impurity model for a range of the discretization parameter A. 



fects in NRG calculations is the "self-energy trick" 6 . It con- 
sists of numerically computing the self-energy as the ratio 

E^w) = UF a (u;)/G a (Lj) where F a (u) = ((n_ CT d CT ; 4)) u 
and G a (uj) = {{d a \ and then computing an improved 
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Green's function as 



VI. DISCRETIZATION ARTEFACTS 



G| T mprovcd (w) 



u> - e- £(w) + A(u>) 



(20) 



An additional merit of this technique is that it leads to a par- 
tial cancellation of the oscillatory features in G a and F a , 
giving a smooth self-energy S CT . In Fig. [3] we compare raw 
and self-energy-improved spectral functions for the symmet- 
ric and asymmetric Anderson model. We first note that the 
change of the spectral function upon using the self-energy 
trick is rather small, unlike in the case of log-Gaussian broad- 
ening with large b where the self-energy trick leads to a siz- 
able improvement and reduction of overbroadening. Results 
for the symmetric case (Fig. |3^) show that while the Friedel 
sum rule is satisfied to better accuracy, the self-energy trick 
leads to slightly broken particle-hole symmetry in the final re- 
sult, which is not desirable. On the other hand, in the general 
asymmetric case the self-energy trick cures problems associ- 
ated with different limiting behavior of A{uS) for uj — ► + and 
u> — > CP, respectively (see Fig.[3]3, inset with the close-up on 
the Kondo resonance). 
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Figure 3: (Color online) Spectral function of a) symmetric and b) 
asymmetric Anderson impurity model: comparison of raw spectral 
function with that obtained using the self-energy trick. 



A. Types of artefacts 

Closer inspection of the computed high-resolution spectral 
functions reveals the presence of artefacts which cannot be 
entirely eliminated by increasing N z or reducing A. These are 
thus genuine intrinsic NRG discretization artefacts. 

As a first example, we plot in Fig. [4^ the spectral func- 
tion A(u>) of the non-interacting impurity model in the high- 
energy range near the band edge, i.e. the tail of the Lorentzian 
spectral peak. We see a pronounced artefact which shifts to- 
ward the band-edge as A is decreased. If the exact solution is 
subtracted from the artefact, we find that there is some can- 
cellation of positive and negative differences, but there is nev- 
ertheless a positive net (integrated) difference; this is the ori- 
gin of the previously mentioned missing spectral weight in 
the spectral peak of the resonant-level model. In the inset we 
show the spectral function Af (u>) of the first site of the Wil- 
son chain /q. For a flat band, p — 1/2D, this function should 
likewise be flat, except for the features that are mirrored from 
the impurity spectral function A(u). The NRG discretization, 
however, introduces additional artefact structure for energies 
near the band-edge. 

To study this problem more closely, we compute A f Q (oj) 
for a system without the impurity, Fig. @J}. In addition to the 
very pronounced band-edge artefact, there are also discernible 
additional artefacts at lower energies. The ratio of energies 
for two consecutive artefacts is A, as expected. The artefact 
peaks presumably exist down to lowest energy scales, but their 
amplitudes decrease rapidly and eventually the peaks can no 
longer be resolved since they are masked by the residual os- 
cillations in the calculated spectral functions. Curiously, the 
average value of Af (u>) in the low-energy region seems to 
have a minimum for A « 1.8. Furthermore, for this value of 
A, the artefacts appear to be the largest. This is in agreement 
with the results for spectral functions presented above. These 
artefacts can, however, be strongly reduced finding proper pa- 
rameter p of the spectrum patching procedure. 

In Fig. UJ; we plot the spectral density Af (uj) for different 
values of p. If p is too small, we obtain very pronounced 
discretization artefacts. If p is too large, the spectral density 
is underestimated. The optimal value of p is around 2, but 
it depends on the energy cutoff in the truncation; we work 
with cutoff E cuto fi = IOwn, thus for p = 2 and A = 2, we 
have pA 2 u>N = 8ojn < E cuto s- We remark that the large 
artefacts near the band-edge are not related to the patching 
procedure (see also below), although the form of the artefacts 
does depend somewhat on the value of p. 

We can formulate the following recipe for choosing appro- 
priate NRG parameters: 

1. fix A; 

2. increase truncation cutoff until NRG results no longer 
change significantly; 

3. tune rj and N z to suppress overbroadening of spectral 
functions; 
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Figure 4: (Color online) a) High-energy artefacts in the spectral function of the resonant-level model. Inset: the spectral function on the first 
site of the Wilson chain, fo. b) Spectral function on the first site of the free Wilson chain, A/ , for different values of the discretization 
parameter A. c) Spectral function A f computed for different values of the spectral patching parameter p. d) Spectral function A / in the case 
of semi-elliptic DOS, p(e) = Poa/1 — (e/D) 2 . e) Spectral function Af in the case of semi-elliptic DOS, p(e) = po\A ~ (2e/-£>) 2 with 
support [—0.5 : 0.5]D. f) Comparison of spectral function Af computed with the complete-Fock-space NRG approach and the conventional 
density-matrix NRG approach. 



4. tune p for good reproduction of the band spectral func- 
tion A fo (u). 

If necessary, steps 2-4 may be reiterated. To be specific, for 
A = 2, E cuto3 = l0u N , T) = 0.01, and N z = 64, we find 
that Af a (u>) is closest to 1/2D at low energies forp = 2.1. A 



caveat is in order: tuning p for good reproduction of Af (u>) 
does not necessarily imply that the same value of p will be 
optimal for the full problem (with the impurity coupled to the 
bath). Nevertheless, such p is most likely a good choice. 

For applications of the NRG as an impurity solver in 
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DMFT, it is important to reproduce an arbitrary conduction- 
band DOS as accurately as possible. As a simple test, in 
Fig. HJi we consider the case of the cosine band dispersion, 
ej; = D cos k, which has a semi-elliptic DOS, 



(21) 



We again find sizable artefacts near band edges at approxi- 
mately the same positions as in the case of a flat band. One 
might expect that using a DOS with a limited support (such 
that it excludes the strong artefacts at w 0.7 D) would resolve 
the issue. Alas, that is not the case. The artefacts simply ap- 
pear at rescaled positions, as is shown in the example of a 
semi-elliptic DOS with support [—0.5 : 0.5], Fig. Any 
abrupt change in the density of states (any sharp feature, in 
fact) is thus expected to lead to anomalies at low energies. 

Spectra calculated using the complete-Fock-space (CFS) 
approach 7 ' 8 also show artefacts, although there are differences 
in the details, see Fig.Hf. There are several advantages to the 
CFS approach: the normalization is satisfied exactly within 
numerical accuracy and there is no ambiguity in the choice of 
the energy range where the spectrum is computed at each it- 
eration (no parameter p). The conventional approach is, how- 
ever, significantly faster since the eigenvectors and matrix el- 
ements need to be computed only in the retained part of the 
Hilbert space in each NRG iteration. In addition, in CFS the 
delta-peak energies are given as a difference between an en- 
ergy of a kept state and an energy of a discarded state; the 
latter is located at the upper end of the shell excitation spec- 
trum, thus it is affected by the accumulated truncation errors 
from all previous NRG iterations. For this reason, the spectra 
calculated using the traditional approach with patching satisfy 
higher-moment sum rules to higher accuracy (in the permil 
range as opposed to the percent range) even though they break 
the normalization sum rule. 



B. Origin of the band-edge artefacts 

In the case of a flat band, p(uj) = const., the origin of 
the main artefact near the band edge is easy to understand. 
Following Ref.Hjl we write the density of states on site /o as 



2D\d£?/dz\ 



(22) 



where ef define the discretization mesh, 



DA 2 - j - 



(j = 2,3,...), 



and £ f are defined as 



3 !j. de/e In (e|/e| +1 ) 



with 



J, which gives 



£f = D 



1 - A" z 
zlnA ' 
1- A- 1 



(25) 



£ i= D i i 
3 In A 



A 2 -^, (.7 = 2,3, 



For given argument u), the parameters z and j in the right hand 
side of Eq. (|22] > are determined by the relation £ = u> which 
has a unique solution. (To simplify the notation and discus- 
sion, we assumed particle-hole symmetry of the conduction 
band and we consider ui > only. All features at positive 
energies are then simply mirrored to negative frequencies.) It 
can be easily shown that for j = 2, 3, . . ., i.e. for 



uj e 



l-A- 
lnA 



1-A- 



lnA 



we indeed have 



A /o (u) = l/2D. 



(26) 



(27) 



This is not the case, however, for j = 1, i.e. for u within 
(1 — A -1 )/ In A from the band edges. We obtain, instead, 



(1 + 0w) 



In A 



(28) 



with /? = W \—e~^ u /J\, where W is the Lambert W- 
function. In Fig. [5] we plot three spectral functions: 1) an- 
alytically calculated spectral function, Aw, 2) the spectral 
function numerically calculated by exact diagonalisations of 
the single-electron Hamiltonians obtained after discretization, 
A^, and 3) the spectral function calculated directly using 
NRG, -4/^ RG) - Compared to the analytical result, AfJ, the 
function A ^ features artefacts due to finite N z and broaden- 
ine, while At' BG ' in addition shows truncation errors. The 

° Jo 

band-edge artefact is thus not some unexpected numerical 
artefact, but it is the direct result of a particular choice of the 
discretization scheme. It arises from a different behavior of 
£1 as a function of z as compared to other £J . This, in turn, 
is due to the presence of the band-edge, which sets the upper 
boundary in the integrals in Eq. (f24-b . 

For arbitrary density of states we introduce weight func- 
tions for different discretization intervals^ 



PjO 



P{e) 



1/2 



In p(u)du; 



(29) 



(23) so that the operator f Q takes the following form 



1/2 



(30) 



(24) where aj m are conduction-band operators for the m-th mode 
(combination of states) in the j-th discretization interval; only 
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Figure 5: Comparison of spectral functions Aj a . 



m — modes are retained in the NRG. The spectral function 
on the first site of the Wilson chain is then given as 



Jj p{uj)duj 
AM = \d£*/dz\ ' 



(31) 



where z and j are again determined by the relation £? = u>. 
In order to achieve decoupling of higher modes in each dis- 
cretization interval, Campo and Oliveira proposed to calculate 
coefficients £ J ap 2 ' 



'* f IjP (e)/ede- 



(32) 



In the most commonly used conventional discretization 
scheme 32 , the coefficients are given, instead, as 



f IjP (e)ede 

L p( £ ) de ' 



(33) 



It is easy to verify that £| calculated in either way do not 
satisfy the equation Af (uj) = p(u>) and that strong artefacts 
appear near sharp features in the density of states. As an ex- 
ample, we compare in Fig. |6]the cosine band DOS with Af a 
computed with both discretization schemes. Both show sig- 
nificant band-edge artefacts (see also Fig. |4}i). In the con- 
ventional scheme, the spectral function A/ in addition sys- 
tematically underestimates p(ui) at lower energy scales by the 
well-known factor of 



A A 



lnAl + A" 
2 1-A- 



(34) 



which is taken into account in practical NRG calculations in 
an ad-hoc manner by multiplying the impurity hybridisation 
(or exchange constant) by this same value. 



VII. OVERCOMING THE DISCRETIZATION ARTEFACTS 

We have demonstrated that the origin of the discretiza- 
tion artefacts is in the z-dependence of the coefficients 
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Figure 6: Analytically computed spectral function Af for semi- 
elliptical DOS in the conventional and Campo-Oliveira discretization 
scheme compared with the exact DOS. 



These coefficients, in turn, are determined by the discretiza- 
tion points €j, the choice of the basis states of the discretized 
conduction band (in particular the weight functions cf>jo) and 
the recipe for the calculation of coefficients 2 -^. Keeping 
the same set of the discretization points and zero-mode ba- 
sis states, we may decide to define £* in a more appropriate 
way, i.e. such that all coefficients satisfy the condition 



\d£*/dz\ 



p{u). 



(35) 



Details about solving this equation are given in the Appendix 
B. Well-behaved solution may be found for arbitrary DOS 
function p(co) and the asymptotic (large j) behavior of £* is 
the same as in the Campo-Oliveira discretization scheme. 

We note that this modification of the discretization proce- 
dure in no way makes NRG an exact method, even though 
we expect much better reproduction of the conduction band 
DOS. In the spirit of the original NRG procedure, we still rely 
on the assumption that discarding higher-mode states in each 
discretization interval is a good approximation which can be 
systematically improved by reducing A toward 1 . In particu- 
lar, discretization-related artefacts are still possible and we in- 
deed find them, as detailed in the following. The improvement 
consists in significantly reducing the severity of the artefacts. 

Solving Eq. ( [35] ) in the case of a flat band, only £f is mod- 
ified, while £* for j > 2 remain the same. We obtain 



£1 



1- A" 
lnA 



l-z. 



(36) 



As z is swept from to 1, this quantity takes values over the 
same interval as the Campo-Oliveira expression for £f . This 
is important, since £* must cover the whole energy range. 
In Fig. [7] we compare the spectral function Af (u>) computed 
with original and modified discretization approach. The im- 
provement is, as expected, significant. The spectral function 
overshoots slightly (by less than two percent) as the band-edge 
is approached and it decays to zero on the scale set by the 
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broadening parameter 77. A closer look reveals small resid- 
ual artefacts positioned at energies £j =1 , J = 1,2,.. ., which 
take the form of asymmetric dips. Their weights rapidly de- 
creases with increasing j; in the worst case, for j = 1, the 
dip amplitude is less than one permil of the background 1/2D 
weight. There are further artefacts between the £j =1 dips, but 
their amplitudes are even smaller than those of the main arte- 
facts. At low energies, Af {u>) converges to 0.50025, which 
can be tuned exactly to 1/2 by further tuning of the patching 
parameter p. 
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Figure 7: Comparison between the Campo-Oliveira (old) and the im- 
proved discretization (new) approaches. 

In Fig.[8]we compare the spectral functions of the resonant- 
level model obtained using both discretization schemes. We 
find that the new discretization scheme strongly suppresses 
the artefact peak structure and correctly reproduces the behav- 
ior at the very edge of the conduction band (within the limits 
imposed by the broadening procedure). We also see that the 
flanks of the spectral peak agree better with the exact solution. 
On the other hand, we see that an artefact appears at the very 
top of the resonance. This artefact is directly connected with 
the discretization itself and does not depend, for example, on 
the truncation or patching; the situation improves, however, 
with decreasing A (see Fig. Q~2] in Subsection B below). We 
point out that the artefact is not located at any £f , thus it 
is not related to the residual artefacts found in Af (u>) of the 
decoupled band. It should rather be interpreted as a finite-size 
effect due to representation of the continuum by a finite chain; 
the ^-averaging cannot entirely eliminate such effects. In spite 
of the artefact, we may conclude that the overall reproduction 
of the spectral weight distribution is considerably improved. It 
may also be noted that we present here the most difficult case: 
a very broad resonance near the band edge. Such situation is 
rather unusual for impurity problems; for narrow resonances 
and for peak energies closer to the Fermi level the double- 
peak artefact is quickly reduced. Broad spectral distributions 
are, however, typical for DMFT applications, where residual 
artefacts may become more problematic. 

In Tables [III] and [IV] we show the moments for the non- 
interacting model and for the asymmetric Anderson model. 
They are to be compared with the corresponding Tables [I] and 
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Figure 8: Spectral function of the resonant-level model: comparison 
between the exact analytical result and two different discretization 
approaches. 



HI1 The agreement between (jf-"' and in the new scheme 
is below one permil for all moments, as in the old one. In 
the resonant-level model, the agreement of calculated /14 now 
agrees with the exact value within one permil (while in the 
Campo-Oliveira scheme we found a discrepancy of 7%). In 
the Anderson model, we also observe a change in ^4 of the 
same order, suggesting a similar degree of improvement. 



A. Arbitrary density of states 

In Fig.[9]we demonstrate on the example of the semi-elliptic 
DOS that the proposed discretization approach can also be ap- 
plied for an arbitrary density of states. In this case, all £? are 
modified and they need to be numerically calculated using the 
technique described in the Appendix B. As in the case of flat 
band, some small discrepancies between Af (u>) and p(co) are 
found at the very edge of the band. The over-all agreement is, 
however, significantly improved on all energy scales. 




Figure 9: Spectral function Af in the case of semi-elliptic DOS 
computed using the new discretization scheme. 

We test the method on the case of a symmetric Anderson 
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Moment 



Exact, p,\ 



(e) 



Static, (i\ 



Dynamic (delta peaks), (i\ 



(d) 



Dynamic (broadened), fi. 



(b) 



(10 
(11 
(12 
(13 
(14 



1 

-0.050000 
0.0056831 
-0.00044331 
0.00110129 



-0.050000 
0.0056831 
-0.00044331 
0.00110120 



0.999979 

-0.049998 

0.0056876 

-0.00044376 

0.00110158 



0.999964 

-0.049997 

0.0056704 

-0.00044217 

0.0010842 



Table III: Moments for the non-interacting impurity model with parameters e/D — —0.05 and T/D = 0.005. Improved discretization 
scheme, A = 2, r? = 0.015, N z = 32, p = 2. 



Moment Static, /i' s ' Dynamic (delta peaks), ^i' d ' Dynamic (broadened), fi\ 

/hi 1.000302 1.000287 

Hi -0.0123213 -0.0123193 -0.0123187 

fi 2 0.00455274 0.00455661 0.0045386 

fi 3 -0.000138138 -0.000138191 -0.000137434 

fi4 0.00109664 0.00109694 0.00107882 



Table IV: Moments for the asymmetric Anderson model with parameters U / D = 0.07, e/D — —0.05, T/D = 0.005. Improved discretization 
scheme, A = 2, r\ = 0.015, N z = 32, p = 2. 



model with semi-elliptic DOS. In Fig. \TU\ we plot spectral 
functions for rather large T — 0.1D for increasing values of 
U. While for small U the functions are rather smooth, we 
observe more pronounced residual artefacts for large values 
of U, as the charge-transfer (Hubbard) peaks approach the 
band edge (see, for example, the U/D =1.5 case). Neverthe- 
less, the results are significantly more physically sensible than 
those obtained using conventional broadening and discretiza- 
tion techniques. For U > 2D, the Hubbard peaks are located 
outside the conduction band. They become narrower and they 
have strongly asymmetric shape 33 ; in fact, in some parameter 
ranges they have a two-peak structure. We also note that the 
impurity parameters used here are comparable to those that 
typically arise in effective models in DMFT (see also Sec.lIXIi. 
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Figure 10: Spectral function of the symmetric Anderson model with 
semi-elliptic DOS computed using the new discretization scheme. 



B. Convergence with A 

The new discretization scheme vastly improves the conver- 
gence to the A — * 1 limit. We demonstrate this in Fig. QT| 
by comparing the calculated level occupancy in the resonant- 
level model with the exact value as a function of A. With 
the new approach, we obtain very accurate results even with 
very large discretization parameter (four digits of accuracy at 
A = 8). In other approaches, not only is the convergence 
to the continuum limit slower, but extrapolating the numeric 
results in the range A > 1.5 to the A — > 1 limit leads to 
a systematic error; presumably the assumption of quadratic 
(or polynomial) A-dependence no longer holds for smaller 
A. With the improved discretization approach one can com- 
pute expectation values of various operators reliably even at 
very large A: this is quite important for numerically demand- 
ing multi-orbital/multi-channel quantum impurity problems. 
Similar improvements also hold for calculations of thermo- 
dynamic quantities (such as the impurity contribution to the 
magnetic susceptibility and entropy). 

We have seen previously that residual artefacts are observed 
in spectral functions. In Fig. [12] we report how these residual 
artefacts are reduced as A is reduced. For sufficiently small A, 
the artefact appearing as double peak structure is eliminated. 
Furthermore, we see that the artefacts shift as a function of 
A. This implies that some additional improvement could be 
obtained by performing the calculation for several different 
values of A and averaging the resulting spectral functions. 

In the sense that the new discretization scheme gives the 
best possible representation of the conduction band DOS by 
the Wilson chain (after the z-averaging), this technique pro- 
vides the best results that one can achieve by representing each 
discretization interval by a single level. A possible system- 
atic improvement would consist in including more than one 



12 



1.516 



1.514 



1.512 



1.51 



->^- V ^J*-»--*- 



U=0, e/D=-0.05, I7D=-0.05 
E {f =10cn 

cutoff N 



• new 

■ Campo-Oliveira 

♦ conventional 
^ exact result 



1.5 



2.5 
A 



3.5 
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Figure 12: Convergence of the spectral function of the resonant-level 
model to the exact result with decreasing A in the case of the new 
discretization scheme. 



mode for low j (where the band DOS still varies strongly as 
a function of energy) and performing the NRG in the star ba- 
sis, perhaps using Lanczos exact diagonalisation procedure to 
diagonalize the cluster at each NRG iteration. 



C. Spectral features outside the conduction band 

We tested how accurately the NRG reproduces spectral fea- 
tures at energies outside the energy band (i.e. outside the 
[— D : D] interval in the case of a flat conduction band) for 
the example of the resonant-level model. In this model, for 
e < —D, there is a <5-peak at the energy cjq given by 



loq - e + ReA(w ) = 0, 



with weight 



1 + 



1 

' dRcA(u) 



(37) 



(38) 



while the spectrum in the [— D : D] range is described by 
Eq. ( fT8l ). We compare the calculated spectrum with the ex- 
pected results in Fig. Qj] The <5-peak takes the form of the 
broadening kernel, Eq. (O, and we can accurately extract its 
position, height and width by fitting to an exponential func- 
tion Aexp[— (u> — luq) 2 /2a 2 ]. We find that the position and 
the (integrated) weight of the peak are reproduced within ap- 
proximately four digits of precision. Furthermore, we find 
that 



a/ajQ = 0.01009, 



(39) 



which is to be compared with the broadening factor r\ = 0.01. 
We conclude that within one percent accuracy, there is no 
other source of broadening than the explicit spectral function 
broadening by the Gaussian broadening kernel. The agree- 
ment of the calculated spectral function within the conduction 
band, i.e. in the [— D : D] interval, Fig. [T3~b. with the exact 
result is also very satisfactory. 
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Figure 13: Spectral function of the resonant-level model in the case 
where the resonance is outside the conduction band. 



VIII. HIGH-RESOLUTION SPECTRAL FUNCTIONS 

We present two examples of high-resolution calculations 
unmasking interesting details. We first study the emergence 
of the Kondo resonance as the electron-electron repulsion U 
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is increased in the Anderson impurity model. We then con- 
sider the Anderson-Holstein model to show that the phononic 
side peaks can be well resolved. 

A. The emergence and the shape of the Kondo resonance 

In Fig.[T4]we show spectral functions for the Anderson im- 
purity model for a range of values of the electron-electron 
interaction U, from the non-interacting case to the symmet- 
ric situation U = — 2e (Fig. fT4h ) and then to the large-[/ 
limit (Fig. [T4"b). Since these results are hardly affected by 
overbroadening, we can accurately follow the evolution of the 
spectral peak, its location as well as its height and width. In 
the non-interacting limit, the peak height is 1 /irT, its width is 
w r and it is centered at u> « e. As U increases, the peak 
position shifts linearly with U (Hartree shift), while its height 
decreases. The remaining spectral weight is located in the 
emerging lower charge-transfer spectral peak (i.e. the lower 
"Hubbard band"); this peak is initially located below e, but it 
shifts to rj e as we approach the particle-hole symmetric situ- 
ation. The width of the charge-transfer peaks is roughly twice 
(2T) the width of the original non-interacting peak (F). As we 
increase U further, Fig. [T4b . the lower charge-transfer peak 
shifts only weakly as a function of U, while the upper charge- 
transfer peak shifts as e + U; in the range of finite U shown, 
its height decreases only slightly and the width remains nearly 
constant. At the same time, the width of the Kondo reso- 
nance is significantly reduced, but we find that it remains al- 
most pinned at the Fermi level (at U — oo, for example, the 
half-width at half-maximum of the Kondo peak is 1.2 10 _8 -D, 
while the shift of the maximum is only 3.6 10~ 10 Z?, i.e. 3 
percent in the units of HWHM). This is in agreement with 
the Fermi liquid theory, but in disagreement with the results 
from methods based on the large- N expansion, such as the 
non-crossing approximation, which overestimate the shift of 
the resonance, in particular for N = 2. It also implies that 
the Kondo temperature should better not be defined as the dis- 
placement of the Kondo resonance from the Fermi level, as it 
is sometimes done. 

In Fig. [T5b we plot a close-up on the Kondo resonance in 
the symmetric case, e + U/2 = 0. As expected, the peak 
shape deviates significantly from a Lorentzian shap o 34 i 3S i 36 ' 37 . 
In fact, true agreement is only found in the asymptotic uj — > 
region, where both the Lorentzian curve and the spectral func- 
tion have quadratic frequency dependence. In the latter case, 
this is mandated by the Fermi-liquid behavior at low energy 
scales. 

The relation between the width of the Kondo resonance and 
the Kondo temperature Tk (times k B ) is of considerable ex- 
perimental interest, in particular for tunneling spectroscopy. 
In the symmetric case, we find for the ratio between the half- 
width at half-maximum and the Kondo temperature (Wilson's 
definition): 

A-hwhm /Tk,w ~ 3.7. (40) 

The Kondo temperature Tk w is defined as Ximp (T — 
0) = (gti B ) 2 (W/4n)l/k B T K , w « (g fi B ) 2 0. 103 /k B T K ,w, 
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Figure 14: (Color online) Spectral functions of the Anderson model 
for increasing U. 




Figure 15: (Color online) Close-up on the Kondo resonance of a) 
symmetric and b) asymmetric Anderson impurity model and a fit to 
a Lorentzian (red curve) in the Fermi liquid regime for uj <C Tk ■ 
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where W = e c+1 / 4 / ^/n is Wilson number, and extracted 
from the NRG results for the magnetic susceptibility using 
the prescription k B T K . w Ximp(T K . w ) / '(ffMs) 2 = 0.07. The 
same value of 3.7 is also obtained when log-Gaussian broad- 
ening is used with a small value of the parameter b and with 
suitable z-averaging. This ratio is lower than some other val- 
ues reported in the literature^. 

In Fig. [T5b we plot the Kondo resonance in the asymmet- 
ric case. We find that the ratio between the half-width and 
the Kondo temperature is now Ahwhm/^V.w = 4.6. Even 
though we are still deep in the Kondo regime (phase shift is 
5 w 0.477r), the Kondo peak has developed a significant asym- 
metry in its shape. These line-shape effects are important in 
the interpretation of experimental results. Due to uncertainties 
in the ratio Ahwhm/TV.w^ the expected systematic error in 
determining Tk from the Kondo peak width is estimated to 
be several 10 percents. This implies that comparisons of Tk 
of different adsorbate/surface systems determined in this way 
are rather meaningless unless the differences are of the order 
of a factor 2 or more. 



B. The phononic side-peaks in the Anderson-Holstein model 

We consider the Anderson-Holstein model with coupling of 
a local Einstein phonon mode to charge fluctuations: 

-Himp = en + U n^ni + g(n — l)(a^ + a) + lo^o) a. (41) 

Here a is the bosonic phonon operator, ujq is the oscillator fre- 
quency and g the coupling between the impurity charge and 
the oscillator displacement. This model was studied intensely 
using a variety of techniques, including NR G 39 ' 40 ' 41 ' 42 ' 43 . 
Its applications range from the problem of small polaron 
and bipolaron formation, electron-phonon coupling in heavy 
fermions and valence fluctuation systems, to describing the 
electron transport through deformable molecules. 

The effect of the electron-phonon coupling is to reduce 
the effective electron-electron interaction and shift the level 

39 44 

energ y^ ' : 



g 2 

e e ff = e H . 

UJo 



(42) 



In addition, the effective hybridisation becomes phonon- 
dependent, since the phonon cloud can be created or absorbed 
when the impurity occupancy changes 39 . 

It is possible to resolve the phononic side-peaks and the 
transition to the charge Kondo regime, Fig. [16] For small 
coupling g, we see the gradual emergence of the phononic 
side-peaks at energies e ff + U c g + nuio, n = 1,2,3,.... 
In addition to these peaks, we see that the charge transfer 
peak at e c ff + U c g itself has internal structure; as g increases, 
part of the spectral weight is transferred from this peak to 
higher energies in the form of a smaller peak which eventually 
merges with the first phononic side-peak at e c g + U c s + o->o- 
The transition from spin to charge Kondo regime occurs at 



g/D w 0.0445, when U e e ss 0. At the transition, the charge 
transfer peak merges with what used to be the Kondo res- 
onance to give a single broad resonance whose width is no 
longer set by the energy scale of the Kondo effect, but rather 
by some renormalized spectral width r c ff . 
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Figure 16: Spectral functions for the Anderson-Holstein model in the 
particle-hole symmetric case, S = e + U/2 — 0. 



IX. NRG AS A HIGH-RESOLUTION IMPURITY SOLVER 
FOR DMFT 

The most severe shortcoming of the NRG (using log- 
Gaussian broadening with large b and traditional discretiza- 
tion schemes) in its role as an impurity solver in DMFT was 
the reduced energy resolution at finite excitation energies. 
This not only affects the self-consistent calculation by intro- 
ducing systematic errors, but sometimes features in spectral 
functions at high energies (for example kinks in the excitation 
dispersions) are themselves of interest. We demonstrate the 
applicability of the new approach on the simplest example of 
the Hubbard model. The case of hypercubic lattice is consid- 
ered in Fig. [T7] where we plot the local density of states for 
a range of the repulsion parameter U as the metal-insulator 
transition is approached. Compared to the results computed 
using the conventional NRG approach, the high-energy fea- 
tures (Hubbard bands) are sharper. Furthermore, the conven- 
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tional approach underestimates the reduction of the density 
of states ("pseudo-gap") between the Hubbard bands and the 
quasiparticle peak. We also observe that the Hubbard bands 
have inner structure. We find a notable peak at the inner edges 
of the Hubbard band; the existence of some spectral features 
at the band edges had been suggested already in the early iter- 
ative perturbation theory, non-crossing approximation, quan- 
tum Monte Carlo and NRG DMFT results for the Hubbard 
model and the existence of a sharp peak was demonstrated 
in more recent high-resolution dynamic density-matrix renor- 
malization (D-DMRG) calculations^^. 
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Figure 17: Local spectral functions of the Hubbard model on the 
hypercubic lattice. 

On the Bethe lattice, Fig.[l8j the Hubbard bands are sharper 
due to the finite support of the lattice density of states and the 
inner Hubbard band edge peaks are sharper. There are further- 
more less pronounced spectral features at integer multiples of 
the energy of the inner Hubbard band edge; they are most vis- 
ible in the U /W = 1.4 results. We also calculated the local 
two-particle Green functions at the end of the DMFT cycle to 
try to obtain some insight whether these additional structures 
are possibly related to certain two-particle excitations. How- 
ever, no clear evidence was found for such statement. Thus, at 
present, we cannot give a satisfactory physical explanation of 
these additional structures. In any case, they motivate further 
high-resolution studies of both the single-impurity Anderson 
model and the Hubbard model in DMFT, concentrating on the 
regime with vanishing Kondo resonance. 



X. CONCLUSION 

We presented spectral function calculations which indicate 
that the numerical renormalization group method allows to 
compute more accurate results than it is generally believed. 
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Figure 18: Local spectral functions of the Hubbard model on the 
Bethe lattice. 



We have shown that overbroadening effects can be in large 
part removed by using a sufficiently narrow Gaussian broad- 
ening kernel. Furthermore, we have shown that there is sur- 
prisingly little variation as A is decreased (disregarding the 
artefact shifts), thus there is no inherent overbroadening due to 
the discretization of the conduction band. At best one can say 
that as A is increased, more z values need to be used in the in- 
terleaved method to obtain smooth spectral functions. It must 
be emphasized that sweeping over z is an embarrassingly par- 
allel problem, i.e. essentially no overhead is associated with 
splitting the problem into a large number of parallel tasks. 

As the continuum limit is approached (A — > 1), the dis- 
cretization artefacts in the spectral function calculated using 
the traditional schemes shift out toward the band-edge, but 
in the range of A that can be used in practical calculations, 
the artefacts are always present. The use of the logarithmic 
discretization is commonly justified by the rapid convergence 
of calculated quantities to the continuum limit; while static 
properties indeed converge rapidly, this is not the case with 
dynamic properties. The presence of artefacts therefore has 
several implications for NRG calculation. First of all, in the 
traditional approach it cannot be claimed that a calculation is 
performed for a given density of states p(uj), but rather for 
a band with a density of states given by Af a (oj) in the prob- 
lem with decoupled impurity. The presence of the structure 
in the spectral function Af (u>) then forcibly leads to what 
is perceived as "artefacts" in the impurity spectral function 
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A(u>). Artefacts have important implications for the appli- 
cation of the NRG in DMFT, since these anomalies lead to 
features in the impurity spectral function that are difficult to 
disassociate from real fine structure. A good test to distin- 
guish between artefacts and real spectral features is to per- 
form calculations for several values of A, keeping all other pa- 
rameters constant. Real features will change very little, while 
artefacts will shift and change form significantly. Depending 
on the circumstances (structure of the impurity model, model 
parameters) and the purposes (single-impurity calculation vs. 
self-consistent dynamical mean-field-theory calculation), the 
artefacts are either benign or rather detrimental. 

The proposed new way of calculating the coefficients £*■ 
leads to a sizable improvement in the convergence to the 
A — > 1 limit and to a significant reduction of the discretization 
artefacts. Since the DMFT self-consistency loop couples low- 
energy and high-energy scales, the reduction of the artefacts 
at high energies is a significant improvement which increases 
the reliability of the NRG as an impurity solver. 
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Appendix A: SPECTRAL FUNCTION BROADENING 



compare the results for spectral functions obtained using dif- 
ferent broadening kernels, Fig. [19] We compare simple Gaus- 
sian broadening, conventional log-Gaussian broadening and a 
modified log-Gaussian kernel proposed in Ref.H: 



P(u,E) 



fin w-hE 



y/*a\E\ 



(Al) 



with 7 = a/4. The same A, p and N z were use for all three 
broadening kernels, with b = a = y2rj. 

In the low-energy (uj <C Tjc) range, we find that 
Gaussian broadening overestimates the spectral density, log- 
Gaussian broadening underestimates it, while the modified 
log-Gaussian kernel, Eq. flAU . gives a very good approxima- 
tion to the high-resolution result. All three approaches de- 
scribe quite well the flanks of the Kondo resonance. Gaus- 
sian broadening overestimates the spectral density in the en- 
ergy range between the Kondo resonance and the Hubbard 
peak, while the best results are here obtained by the original 
log-Gaussian broadening. All three broadening approaches 
shift the maximum of the Hubbard peak to lower energies 
to roughly comparable degree. Finally, in the high-energy 
range, log-Gaussian approaches overestimate the spectral den- 
sity more than the simple Gaussian broadening. 

For studying low-energy properties with typical NRG 
broadening parameters, the modified log-Gaussian kernel, 
Eq. (IA It . is the best choice. For high-resolution studies with 
very small broadening, all three broadening techniques be- 
come almost equivalent, but the plain Gaussian kernel has a 
small advantage by being symmetric; the symmetry leads to 
smaller deviations of higher-moment spectral sum rules. 



Appendix B: MODIFIED DISCRETIZATION SCHEME 




Figure 19: Spectral function of the symmetric Anderson impurity 
model: comparison of results obtained using different broadening 
kernels. Reference results are calculated using Gaussian broadening 
with sufficiently narrow kernel so that very little change is obtained 
by further narrowing. 

In practical NRG calculations, z-averaging is performed 
over a smaller number of twist parameters, therefore wider 
broadening functions must be used. It is thus interesting to 



We describe the modified discretization scheme which con- 
sists of solving the ordinary differential equation for £? : 



4 P( e ) de 
\d£ z Jdz\ 



p(ui). 



(Bl) 



As a first step, we introduce continuous indexing as x = j + z 
with parameter x running from 1 to +oo, so that coefficients 
£j and e| become continuous functions of x, i.e. £(x) and 
e(x). We then rewrite Eq. ( IB 1 b as 



d£{x) 
dx 



re(x+l) 
Je(x) 



p(uj)du> 



P [£{x)] 



(B2) 



with the initial condition £ (1) = D. It is helpful to take into 
account the expected asymptotic behavior of £ (x) using the 
following Ansatz: 



£{x) = Df(x)A 2 -*, 
with /(l) = 1/A. The equation to solve is then 

(x) 



df(x) 
dx 



In A f{x) 



A 2 -x p {£{x)} ' 



(B3) 



(B4) 



17 



This equation can be solved numerically; for general p{u>) it is level, we must have /(oo) = (1 — A -1 )/ In A. Checking 
advisable to use arbitrary-precision numerics for this purpose, the convergence to this value is a good test of the integration 
since the equation is stiff. For DOS which is finite at the Fermi procedure. 
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